import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500),
hv.opts.Image(width=500, colorbar=True, cmap='Viridis'))
import numpy as np
import scipy.signal
import scipy.fft
from IPython.display import Audio
6. Espectrograma¶
6.1. ¿Qué ocurre con el espectro de una señal si su frecuencia cambia en el tiempo?¶
Consideremos este ejemplo sencillo de una señal donde la frecuencia cambia abruptamente en un tiempo determinado
Sea por ejemplo \(f_1=440\) Hz y \(f_2 = 220\) Hz. Si la graficamos:
f1, f2, Fs = 440, 220, 44100
t = np.arange(-0.5, 0.5, step=1/Fs)
N = len(t)
s = np.concatenate((np.cos(2.0*np.pi*f1*t[:N//2]),
np.cos(2.0*np.pi*f2*t[N//2:])))
hv.Curve((t, s), 'Tiempo [s]', 'Señal').opts(xlim=(-0.05, 0.05))
y si la escuchamos
Audio(s, rate=Fs)
Si calculamos la FFT de la señal para obtener su espectro tenemos que
f = scipy.fft.rfftfreq(n=N, d=1./Fs)
S = np.absolute(scipy.fft.rfft(s))
hv.Curve((f, S), 'Frecuencia [Hz]', 'Espectro').opts(xlim=(0, 1000))
Es decir que la DFT/FFT nos entrega un “resumen” de todas las frecuencias en la señal
Importante
No es posible diferenciar una señal donde ambas frecuencias ocurren al mismo tiempo de otra donde las frecuencias aparecen en tiempos distintos
6.2. Frecuencia instantanea¶
Hasta ahora hemos estudiando señales donde la frecuencia es constante en el tiempo
Definamos la frecuencia instantanea como la tasa de cambio del ángulo (fase) en función del tiempo
Una señal sinusoidal con una frecuencia que cambia en el tiempo sería entonces
Notemos que si la frecuencia fuera constante, es decir \(f(t) = f\) \(\forall t\), entonces \(\int_0^t f(\tau) d\tau = t f\) y recuperamos \(A\cos(2\pi t f + \phi)\)
A continuación tres ejemplos de señales donde la frecuencia cambia con el tiempo
El chirrido: señal cuya frecuencia cambia entre dos valores
El vibrato: señal que está modulada en frecuencia por otra señal
6.2.1. Chirrido o Chirp¶
Loica (Sturnella loyca). Referencia: http://www.conserva.cl/2009/09/sonidos-de-aves-de-chile-loica.html
Un chirp es una señal cuya frecuencia varía suavemente entre un primer valor \(f_0\) y un segundo valor \(f_1\). Por ejemplo esta variación podría seguir una forma lineal
donde \(t_0\) y \(t_1\) son los tiempos en que la señal oscila a \(f_0\) y \(f_1\), respectivamente. También se puede usar una forma no lineal, por ejemplo cuadrática o exponencial
Los chirp se usan como modelo en aplicaciones asociadas a radar y sonar. También se han usado para modelar el canto de algunas aves con el objetivo de hacer identificación automática
Podemos crear un chirrido sintético con scipy usando
scipy.signal.chirp(t, # Un vector de tiempos
f0, # La frecuencia en el tiempo t=0
t1, # El tiempo en el cual f=f1
f1, # La frecuencia para el tiempo t=t1
method='linear', # Otras opciones disponibles: 'quadratic', 'logarithmic' o 'hyperbolic'
...
)
f0, f1, Fs = 4000, 2000, 44100
t = np.arange(0, 0.5, step=1./Fs);
s = 0.1*scipy.signal.chirp(t, f0=f0, f1=f1, t1=t[-1], method='quadratic')
En este ejemplo la frecuencia cambia cuadraticamente
hv.Curve((t, f0 + (f1 - f0)*(t/t[-1])**2), 'Tiempo [s]', 'Frecuencia [Hz]')
El resultado sonoro se muestra a continuación
Audio(s, rate=Fs, normalize=False)
6.2.2. Frecuencia Modulada (FM)¶
La FM es una tecnología para guardar información en la frecuencia de una onda electromagnética. Es un tipo de codificación que se usa mucho en transmisiones de radio.
La onda electromagnética se llama señal portadora. En radio corresponde a una sinusoide con una frecuencia central en el rango de 88 a los 108 [MHz]
La información se llama señal modulada. En radio corresponde tipicamente a una señal de voz o a una canción, es decir que está en el rango de los 20 [Hz] a 20 [kHz] (rango audible humano)
Una señal en el rango audible puede viajar algunos metros. En cambio, si va codificada en la señal portadora puede viajar cerca de 50 km
El siguiente esquema muestra la operación que realiza una estación de radio que transmite señales
La radio que recibe la señal debe realizar el proceso inverso, es decir decodificar la información a partir de la frecuencia de la señal que está recepcionando
Matemáticamente la señal modulada \(s_m(t)\) modifica la frecuencia central \(f_c\) de la señal portadora como sigue
donde \(K\) es el coeficiente de modulación y \(s(t)\) es la señal que finalmente viaja por el medio
Cada estación de radio transmite su información \(s_m(t)\) usando una frecuencia portadora \(f_c\) distinta para no traslaparse
6.2.3. Vibrato¶
Un vibrato es un efecto musical que consiste en variar periódicamente el tono de una nota.
Un violinista logra este efecto presionando una cuerda y luego moviendo su dedo de forma regular como muestra la siguiente animación (en cámara lenta)
Podemos considerar el vibrato como un caso particular de modulación de frecuencia. Si consideremos sólo tonos puros podríamos definir \(s_m(t) = \cos(2\pi f_m t)\), con lo que nos queda la siguiente señal
De la expresión tenemos que
\(f_c\) es la frecuencia o tono central
\(f_m\) es la velocidad a la que cambia el tono central
\(K/f_m\) es la amplitud del cambio del tono cnetral
Podemos implementar un vibrato usando
A_c, K, f_c, f_m, Fs = 1, 50, 220, 8, 44100
t = np.arange(0, 2, step=1/Fs)
sm = np.cos(2.0*np.pi*f_m*t)
s = A_c*np.cos(2.0*np.pi*f_c*t + (K/f_m)*np.sin(2.0*np.pi*f_m*t))
La frecuencia de la portadora (azul) aumenta con la amplitud de la modulada (roja)
p1 = hv.Curve((t, s), 'Tiempo[s]', 'Señal', label='Portadora').opts(alpha=0.75)
p2 = hv.Curve((t, sm), 'Tiempo[s]', 'Señal', label='Modulada')
(p1 * p2).opts(hv.opts.Curve(xlim=(0, 0.2)))
Audio(s, rate=Fs)
6.3. Representación en tiempo y frecuencia¶
Para estudiar una señal cuya frecuencia cambia en el tiempo debemos estudiar la evolución temporal de su espectro. La herramienta más utilizada para esto se llama espectrograma
El espectrograma es una representación visual de la energía de la señal distribuida tanto en el tiempo y en la frecuencia. Es decir que es una representación bidimensional.
La siguiente imagen muestra un espectrograma de una señal de habla humana, una señal altamente no estacionaria cuya frecuencia puede presentar cambios bruscos
Notar que:
El eje horizontal representa tiempo (segundos)
El eje vertical representa frecuencia (Hz)
Se usa color para representar la intensidad energética
En la imagen se puede apreciar como el contenido energético cambia su distribución de forma notoria en los momentos de respiración (breath) y habla (speech). Muchos algoritmos actuales de reconocimiento de habla (por ejemplo redes neuronales artificiales) operan reconociendo patrones a partir del espectrograma.
Importante
A diferencia del espectro, el espectrograma nos permite estudiar los cambios de energía instantaneos de la señal
6.3.1. ¿Cómo se obtiene el espectrograma?¶
Para calcular el espectrograma se utiliza la short-time Fourier transform (STFT). Para el caso de una señal discreta la STFT se define como
Notemos como la STFT tanto del tiempo (índice m) como de la frecuencia (índice k)
En la práctica la STFT consiste en
multiplicar la señal por una ventana localizada \(w[n-m]\)
calcular la FFT sobre esa ventana
Esto se repite para distintas ventanas como muestra el siguiente diagrama.
En la parte superior de la imagen la linea azul es la señal y las lineas rojas son las ventanas desplazadas. En la parte inferior se muestra que de cada ventana desplazada se obtiene un espectro. Finalmente el espectrograma consiste en juntar los espectros de amplitud obtenidos para cada ventana. Notemos que las ventanas pueden traslaparse para mejorar la resolución temporal del espectrograma
6.3.2. Ventanas y espectrogramas¶
Para calcular el espectrograma debemos seleccionar
un tipo de ventana
un largo de ventana
un traslape de ventana
Como vimos en la lección de “Fuga espectral” el tipo y largo de la ventana tienen influencia en el resultado
Una ventana debería ser suave para evitar agregar lóbulos laterales
Una ventana debe ser lo suficientemente larga para no perder resolución en frecuencia
Una ventana debe ser lo suficientemente corta para no perder resolución temporal
6.4. Espectrograma en Python¶
Podemos usar la función de scipy.signal.spectrogram cuyos parámetros más relevantes son
spectrogram(x, # Señal
fs=1.0, # Frecuencia de muestreo
window=('tukey', 0.25), # Tipo de ventana y parámetros de ventana
nperseg=None, # Largo de la ventana en número de muestras
noverlap=None, # Cantidad de traslape, por defecto es 1/8 del largo de ventana
...
)
Esta función retorna una tupla con
Un arreglo con las frecuencias del espectrograma de largo M
Un arreglo con los tiempos de las ventanas de largo N
Una matriz de MxN con los valores del espectrograma
A continuación veremos dos ejemplos
6.4.1. Espectrograma de un vibrato¶
Implementemos nuevamente el vibrato con frecuencia instantanea
A_c, K, f_c, f_m, Fs = 1, 25, 440, 8, 44100
t = np.arange(0, 1, step=1/Fs)
sm = np.cos(2.0*np.pi*f_m*t)
s = A_c*np.cos(2.0*np.pi*f_c*t + (K/f_m)*np.sin(2.0*np.pi*f_m*t))
window = ('kaiser', 6.)
A continuación se muestran tres espectrogramas con distinto largo de ventana. La linea roja punteada corresponde \(f(t)\) en función de \(t\). En todos los espectrogramas se usó una ventana de Kaiser con \(\beta=6\)
def plot_spectrogram(Nw):
display(f'Ventana de {1000*Nw/Fs:0.4f} [ms]')
freqs, times, Sxx = scipy.signal.spectrogram(s, fs=Fs, nperseg=Nw, noverlap=Nw//1.5, window=window)
formula = hv.Curve((t, f_c + K*sm)).opts(color='r', line_dash='dashed')
spectrogram = hv.Image((times, freqs, Sxx), kdims=['Tiempo [s]', 'Frecuencia [Hz]']).opts(ylim=(300, 600))
return spectrogram * formula
plot_spectrogram(512)
'Ventana de 11.6100 [ms]'
Ventana muy angosta: Resolución temporal superior (eje horizontal) pero gran dispersión en la frecuencia (eje vertical)
plot_spectrogram(8192)
'Ventana de 185.7596 [ms]'
Ventana muy ancha: Resolución frecuencial superior (eje vertical) pero gran dispersión en el tiempo (eje horizontal)
plot_spectrogram(2048)
'Ventana de 46.4399 [ms]'
Compromiso: La ventana de 46 [ms] parece presentar un mejor compromiso entre resolución temporal y frecuencial.
Importante
El mejor compromiso depende de la tasa de cambio temporal y frecuencial de la señal en particular. No existe una receta para escoger el tamaño de ventana. Lo mejor es siempre probar y estudiar los resultados.
6.4.2. Espectrograma de una señal de voz¶
Consideremos la siguiente señal de voz humana y su espectrograma
Utilizaremos la librería librosa para cargar el archivo de audio en memoria
import librosa
data, Fs = librosa.load("../../data/123.ogg")
time = np.arange(0.0, len(data)/Fs, step=1/Fs)
display(Audio(data, rate=Fs))
hv.Curve((time, data), 'Tiempo [s]', 'Señal')
Nw = 2048
freqs, times, Sxx = scipy.signal.spectrogram(data, fs=Fs, window=('kaiser', 6),
nperseg=Nw, noverlap=Nw//1.5)
hv.Image((times, freqs, 10*np.log10(Sxx+1e-10)),
kdims=['Tiempo [s]', 'Frecuencia [Hz]']).opts(ylim=(0, 2000))
Si comparamos con la señal de audio que graficamos antes podemos observar que
Cada vocal tiene un tono y una distribución de energía distintos
A diferencia de una señal sintética la voz humana es rica en armónicos
A diferencia de una señal sintética hay ruido blanco que contamina el espectrograma
Notar que no graficamos la energía si no su logaritmo en base diez. Aplicar log10 es muy usual para visualizar el espectro en señales de habla ya que su energía varía en un rango muy amplio